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SUMMARY 


A unique dynamic delamination buckling and delamination propagation analy- 
sis capability has been developed and incorporated into a finite element com- 
puter program. This capability consists of: (l)a modification of the direct 

time integration solution sequence which provides a new analysis algorithm that 
can be used to predict delamination buckling in a laminate subjected to dynamic 
loading, and (2) a new method of modeling the composite laminate using plate 
bending elements and multipoint constraints. This computer program is used to 
predict both impact induced buckling in composite laminates with initial delam- 
inations and the strain energy release rate due to extension of the delamina- 
tion. It is shown that delaminations near the outer surface of a laminate are 
susceptible to local buckling and buckling-induced delamination propagation 
when the laminate is subjected to transverse Impact loading. The capability 
now exists to predict the time at which the onset of dynamic delamination buck- 
ling occurs, the dynamic buckling mode shape, and the dynamic delamination 
strain energy release rate. 


INTRODUCTION 

Composite laminates are subject to delamination, which causes a loss of 
both stiffness and strength. Delamination is generally induced by static, 
dynamic, or fatigue loading. Delaminated sublaminates are particularly suscep- 
tible to dynamic local buckling when subjected to Impact loading. The predic- 
tion of impact-induced dynamic delamination buckling is necessary for 
evaluating the durability of many composite structures and Is, therefore, the 
topic of this paper. 

The delamination buckling phenomenon has been observed experimentally 
under static and fatigue loading conditions (refs. 1 to 4), and several analyt- 
ical methods have been proposed to model this damage mechanism. One- 
dimensional beam models and fracture mechanics approaches (refs. 5 and 6) have 
been used to gauge the stability of delaminations in compressi vely loaded lami- 
nates. Finite element approaches (refs. 7 and 8) are often used for these 
analysis. 

In the course of earlier research, experimental observations of dynamic 
delamination buckling in transversely impacted laminates were reported (ref. 9) 
based on high-speed photography of the delamination buckling sequence. This 
work motivated the present development of a finite element analysis technique 
to predict the occurrence of impact-induced dynamic buckling in laminates with 
delaminations. It has been shown (refs. 9 and 10) that the mechanism causing 



extension of a delamination depends on the location of that delamination within 
the impacted laminate. A delamination along the midplane of a symmetrical lam- 
inate will initially extend in a predominantly Mode II fashion. Delaminations 
off the midplane (closer to the outer surface) of a laminate, however, are sus- 
ceptible to buckling caused by compressive stress in the delaminated region. 

This local buckling may then Induce a Mode I dominated extension of the 
delamination. In the latter case, the ability to predict the onset of dynamic 
delamination buckling is essential to determine if extension of the delamina- 
tion will occur when the laminate is subjected to a given Impact load. The ob- 
jectives of this paper are: (1) to outline the computational procedure for 

dynamic delamination buckling and delamination propagation, and (2) to present 
typical results of this procedure for a delaminated composite laminate under 
impact loading. 


DYNAMIC BUCKLING ANALYSIS 

To perform the dynamic delamination buckling analysis, the direct time 
integration solution sequence in the finite element program is altered so that 
a linear buckling analysis is performed at each time step. The buckling analy- 
sis requires solutions of the eigenvalue problem: 


CK] + XCK c ] 


{♦} - 0 


<1) 


where 

CK3 is the structural stiffness matrix 
[K 0 ] is the stress stiffness matrix 

The formulation of these matrices for the NASTRAN Quad-4 plate element used 
here is given in reference 11. 

Each scalar eigenvalue satisfying equation (1) physically represents the nondl- 
mensional ratio 


X - 


q A 
P* 


where a is the time-dependent compressive longitudinal stress in the delami- 
nated sublaminate, A Is the cross-sectional area of the sublaminate, and P* 
is the critical compressive load that will cause buckling of the sublaminate. 
Nhen an eigenvalue reaches the critical value of unity <oA = P*), buckling In 
that mode occurs. The eigenvectors {<j>} associated with each eigenvalue are 
the corresponding dynamic buckling mode shapes. Figure 1 is a flowchart of 
the modified solution procedure. This altered finite element solution proce- 
dure was Implemented in Version 65A of MSC/NASTRAN using NASTRAN DMAP alters. 
The complete modification Is shown in the NASTRAN DMAP Alter sequence In 
Appendix A. 
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FINITE ELEMENT MODELING PROCEDURE 

Figure 2 Is a schematic diagram of the beam-like unidirectional [Olios 
graphite/epoxy laminate with an Initial 5.08 cm (2.0, In.)* long delaminatlon 
through the width and located halfway between the beam midplane and outer sur- 
face. The finite element model of this specimen consisted of a uniform mesh 
of 800 four-node isoparametric plate-bending elements (QUAD-4 elements in 
MSC/NASTRAN) arranged as shown schematically In figure 3. Four elements were 
used In the thickness direction, so e*ch finite element represents 5 plies of 
the 20-ply unidirectional laminate. Although a relatively fine, uniform 
finite element mesh was used for simplicity in this example, a somewhat coarser 
mesh Is allowable in regions where high stress gradients are not anticipated. 

To make the most efficient use of computer time, larger elements may be used 
In regions remote from the delaminatlon without sacrificing any accuracy. 

The allowable displacements of adjoining nodal points In the thickness 
direction were constrained, using multipoint constraint equation sets, so as 
to satisfy simple beam bending assumptions away from the delaminated region 
and In each of the sublaminates In the delaminated region. These assumptions 
are: 


(1) Plane sections remain plane 

(2) No strain In the transverse (Z) direction 

The first constraint set restricts neighboring nodal points in the thickness 
direction to deform along a straight line over the cross section (shear defor- 
mation is neglected) while the beam midplane remains unstrained. The second 
constraint set prevents neighboring plies from penetrating Into each other. 
These constraints are relaxed near the delaminated region so that the displace- 
ments of the nodal points defining the delamination are allowed to deform as 
the solution dictates, independent of the deformation of the neighboring 
plies. This allows the delaminated region to separate from the main laminate 
when a local compression occurs In that area, as shown In figure 4. The con- 
straints used here are shown schematically In figure 5. 

Typical unidirectional graphite/epoxy material constants were used In the 
finite element analysis; E] = 134.4 GPa (19.5x10° psi), E 2 * 10.3 GPa 
(1.5x10^ psi). G 12 = 5.0 GPa (0.725x10^ psi), x>i 2 * °13 * u 23 = 0.3, and 
p = 1580 kg/nw (1.49xl0~ 4 lb-s 2 /1n. 4 ). The 1-dlrectlon (longitudinal) Is as- 
sumed to be aligned with the fibers, while the 2-dlrection (transverse) Is per- 
pendicular to the fibers and in the plane of the laminate. The 3-dlrectlon is 
perpendicular to the plane of the laminate. A representative triangular load- 
ing, shown In figure 2, is applied at the midpoint of the laminate to simulate 
transverse Impact. Both ends of the laminate are assumed to be rigidly 
supported. 


RESULTS AND DISCUSSION 

Figures 6 and 7 show the midplane displacement and bending stress re- 
sponse of a CO] i os laminate without an initial delaminatlon. It is apparent 


‘Original measurements were in inches. 


3 


that the boundaries (clamped ends) have no effect on the structural response 
before 250 ps for the specimen geometry studied. Furthermore, If a delamina- 
tlon exists In the region of high compressive flexural stress 5 to 8 cm (2 to 
3 In.) from the Impact point, local buckling of the sublamlnate region may 
occur In the 150 to 250 ps Interval. 

A 5.08 cm (2.0 In.) long, off-midplane delamlnatlon Is now Introduced 
Into the finite element model as shown In figure 4 by removing the constraints 
described between two adjacent plies In the region from x - 5.08 to 
x = 10.16 cm (x * 2 to x - 4 in.) from the Impact point. The calculated 
dynamic displacement of the partially delaminated specimen under the same load- 
ing Is shown In figure 8. Only the displacement of the right half of the lami- 
nate Is shown, In order to more clearly highlight the opening profile of the 
delamlnatlon. Results of the simultaneous dynamic buckling analysis are shown 
In figure 9. The lowest eigenvalue reaches the critical value of 1.0 approxi- 
mately 190 ps after Impact, for the Impact loading condition shown In figure 2. 
These results Indicate that a first-mode buckling Is likely to occur approxi- 
mately 190 ps after the Impact event begins, for the given loading conditions. 
The buckling may then Induce extension of the delamlnatlon In one direction or 
the other. 

Progressive extension of the Initial delamlnatlon Is modeled by removing 
the constraints between neighboring plies at the appropriate nodal points. 

This approach Is used here to perform finite element analyses of Impacted lami- 
nates with successively Increasing delamlnatlon crack lengths. A simple tech- 
nique Is then used to calculate the variation of the dynamic strain energy 
release rate with time during the Interval on which the Impact load acts. 

Using the applied force, F(t), and the displacement, y(t), at the point of 
Impact, the amount of energy stored In the structure at any time t Is 

r y<t) 

W(t) - J 0 F(£) d£ 

This calculation Is performed for several successive Increments In crack 
length, each Increment being equal to one element length In the finite element 
model . 

Defining 


Gt 



W(a„ + Aa> - W(aJ 
o o_ 

AA 


as the total strain energy release rate; 
where 

a = crack length 
A = crack area 

the variation of strain energy release rate with time at each tip of the 
5.08 cm (2.0 in.) long Initial delamlnatlon Is shown In figure 10. This indi- 
cates that, In the absence of flexural wave reflections from the boundaries. 
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the delami nation will initially e x tend / toward the point of impact. Variation 
of the energy release rate as the initial crack length is extended toward the 
point of impact is shown In figure 11. The apparent trend for the magnitude 
of the energy release rate to decrease as the crack length is Increased indi- 
cates that the available strain energy for delamination extension will decrease 
monotonical ly, leading to stable delamination growth and eventual arrest. 


RELEVANCE TO EXPERIMENTAL OBSERVATIONS 

The dynamic delamination buckling analysis procedure presented here was 
used to predict the observed delamination buckling of the cantilevered speci- 
men shown schematically in figure 12, which was tested as described in refer- 
ence 9. The impact force history and material properties used in the analysis 
are from reference 9. Figure 12 snows the calculated variation of the lowest 
eigenvalue with time, and the associated buckling mode shape. The eigenvalue 
behavior indicates a first-mode del ami nation buckling approximately 150 jiS 
from the beginning of Impact. This is in reasonable agreement with the photo- 
graphic data in reference 9, which shows delamination buckling and initial 
crack extension occurring between 2 to 3 frames <125 to 187.5 ps) from the 
beginning of impact. Further experimental Investigation is currently being 
conducted to provide additional verification and will be reported In the near 
future. 


POSSIBLE EXTENSIONS 

The methods described here can be extended in a straightforward manner to 
analyze the problem of dynamic buckling and extension of an arbitrarily shaped 
embedded delamination in a laminated plate. This problem has been analyzed 
previously under static loading conditions (refs. 7 and 12). The analysis pro- 
cedure for the plate will be similar to that used for the beam-like laminate 
geometry used here, although the physical meaning of the buckling eigenvalue 
will be slightly modified due to the more complex stress field. Further appli- 
cation of this crack-extension modeling procedure to analyze the growth of 
both transply cracks and edge cracks in composite laminates is possible. 
Embodying these extensions into the computational procedure will make it possi- 
ble to assess damage tolerance and structural fracture in composites subjected 
to impact loads. 


CONCLUSIONS 

A dynamic delamination buckling analysis procedure for composite laminates 
with delaminations has been developed and Incorporated into a finite element 
computer program. This algorithm enables the modified finite element program 
to predict: (1) the time at which delamination buckling occurs, and (2) the 

dynamic buckling mode shape. Preliminary experimental verification has sup- 
ported the validity of the analysis. 

The procedure used here to model the composite laminate allows progres- 
sive extension of the delamination to be simulated by removing the appropriate 
constraints between displacements at neighboring nodal points. This slightly 
Increases the total number of degrees of freedom to be solved for in the finite 
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element analysis, while keeping the number of elements and the element mesh 
constant as the crack extends. This procedure greatly simplifies the simula- 
tion of progressive crack extension. It was used to calculate the dynamic 
delamlnatlon strain energy release rate for a composite laminate under Impact 
loading. 

Results obtained show that: (1) delamlnatlons near the outer surface of 

a laminate that Is subjected to transverse Impact loading are susceptible to 
dynamic buckling and buckllng-lnduced propagation, and (2) In the absence of 
boundary effects, transverse impact at a point near a delamlnatlon will cause 
the delamlnatlon to propagate toward the point of Impact, with little or no 
extension In the opposite direction. 

Further application of this procedure to the problem of buckling and prop- 
agation of an embedded delamlnatlon In a laminated plate Is possible. In addi- 
tion, similar techniques can be used to model different flaw types such as 
transply cracks and edge cracks. With these modifications It would be possi- 
ble to computationally assess the damage tolerance of composite structures 
with Initial or Induced defects when subjected to Impact loading. 
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Address Module 


Function 


5 

LABEL 

Begin buckling analysis 

9 

MATMOD 

Strip displacement vector {S} from transient solu 
at each time step 

10 

EMG 

Generate element stress stiffening matrices, [k CT ] 

11 

12-22 

EMA 

Assemble global stress stiffening matrix, [K a ] 
Eliminate boundary conditions; reduce [K a ] 

24 

ADD 

Multiply [K ct ] by (-1 .0) 

25 

DPD 

Red eigenvalue extraction data from input 

27 

READ 

Eigenvalue extraction: 


matrix 


solve 


CK3 


XCK a ] 


{<t>} = 0 for eigenvalues and eigenvectors 


29 

OFP 

Format and print eigenvalues 

30 

SDR1 

Recover dependent degrees of freedom for eigenvectors 

31 

SDR2 

Prepare eigenvectors for output 

32 

OFP 

Format and print eigenvectors 

37 

REPT 

Return to address 5 



FIGURE 1. - DYNAMIC BUCKLING ANALYSIS FLOW CHART. 
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FIGURE 2. - IMPACT SPECIMEN GEOMETRY AND LOADING. 


r 2 ELEMENTS 



FIGURE 3. - FINITE ELEMENT MODEL OF IMPACT FIGURE 4. - IMPACT SPECIMEN AND 

SPECIMEN. BUCKLED SUBLAMINATE. 
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POSITION, CM 

FIGURE 8. - DISPLACEMENT PROFILES OF RIGHT HALF OF IMPACT SPECIMEN WITH INITIAL OFF-MIDPLANE 
DELAMINATION DUE TO A CENTRAL IMPACT. 
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FIGURE 9. - CALCULATION OF DYNAMIC BUCKLING 
INSTABILITY POINT. 





FIGURE 10. - STRAIN ENERGY RELEASE 
RATES FOR EXTENSION OF ORIGINAL 
DELAMINATION IN OPPOSITE DIREC- 
TIONS. 
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FIGURE 11. - STRAIN ENERGY RELEASE RATES FOR 
PROGRESSIVE EXTENSION OF ORIGINAL DELAMI- 
NATION TOWARD THE POINT OF IMPACT. 


FIRST BUCKLING MODE AT 
150 msec 

3.3 



TIME. MSEC 

FIGURE 12. - LOWEST BUCKLING EIGENVALUE 
FOR CANTILEVERED LAMINATE WITH AN 
INITIAL DELAMINATION. 
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